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Abstract 

An approach to pairing in finite nuclei at nonzero temperature is proposed, which incorporates 
the effects due to the quasiparticle-number fluctuation (QNF) around Bardeen-Cooper-SchriefFer 
(BCS) mean field and dynamic coupling to quasiparticle-pair vibrations within the self-consistent 
quasiparticle random-phase approximation (SCQRPA). The numerical calculations of pairing gap, 
total energy, and heat capacity were carried out within a doubly folded multilevel model as well 
as realistic nuclei ^^Fe and ^^'^Sn. The results obtained show that, under the effect of QNF, in 
the region of moderate and strong couplings, the sharp transition between the superconducting 
and normal phases is smoothed out, resulting in a thermal pairing gap, which does not collapse 
at the BCS critical temperature, but has a tail, which extends to high temperature. The dynamic 
coupling of quasiparticles to SCQRPA vibrations significantly improves the agreement with the 
results of exact calculations and those obtained within the finite-temperature quantal Monte Carlo 
method for the total energy and heat capacity. It also causes a deviation of the quasiparticle 
occupation numbers from the Fermi-Dirac distributions for free fermions. 
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I. INTRODUCTION 



Pairing phenomenon is a common feature in strongly interacting many-body systems 
ranging from tiny ones such as atomic nuclei to very large ones such as neutron stars. 
Because of its simplicity, the Bardeen-Cooper-Schrieffer (BCS) theory which explains 
the conventional superconductivity, has been widely employed as the first step in nuclear 
structure calculations that include pairing forces. In infinite systems such as low-temperature 
superconductors, the BCS theory offers a correct description of the pairing gap as functions 
of temperature T and pairing-interaction strength G. Here, as T increases, the BCS gap 
decreases from its value A(0) at T = until it collapses at a critical temperature = 
0.567A(0), at which the phase transition between the superconducting phase and normal 
one (SN-phase transition) occurs 0, sl. However, the application of the BCS theory to 
small systems such as atomic nuclei needs to be carried out with a certain care since quantal 
and thermal fluctuations are not negligible in finite systems, especially when the number of 
particles is small. 

The effects of thermal fluctuations on the pairing properties of nuclei have been the sub- 
ject of numerous theoretical studies in the last three decades. In the seventies, by applying 
the macroscopic Landau theory of phase transitions to a uniform model, Moretto has shown 
that thermal fluctuations smooth out the sharp SN phase transition in finite systems [4] . In 
the eighties, this approach was incorporated by Goodman into the Hartree-Fock-Bogoliubov 
(HFB) theory at finite temperature [sj] to account for the effect of thermal fluctuations Q]. 
Theoretical studies within the static-path approximation (SPA) carried out in the nineties 
also came to the non- vanishing pairing correlations at finite temperature 7|, which are qual- 
itatively similar to the predictions by Landau theory of phase transitions. The shell-model 
and Monte-Carlo shell-model calculations also show that pairing does not abruptly 

vanish at T^,, but still survives at T > T^. For rotating systems, Frauendorf and collabo- 
rators have recently shown a phenomenon of pairing induced by temperature |lfl], which 
reflects strong fluctuations of the order parameter in very small systems with a fixed num- 
ber of particles. The recent microscopic approach to thermal pairing, called modified-HFB 
(MHFB) theory ll|, includes the quasiparticle-number fluctuation (QNF) in the modified 



single-particle density matrix and particle-pairin 
action G is the modified BCS (MBCS) theory |jj 



;ensor. Its limit of constant pairing inter- 
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15| . The MBCS theory predicts a 
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pairing gap, which does not collapse at Tc, but monotonously decreases with increasing T, 
in qualitative agreement with the predictions by the Landau theory of phase transitions and 
SPA. This feature also agrees with the results obtained by averaging the exact eigenvalues 
of the pairing problem over the canonical ensemble with a temperature-dependent partition 
function . The recent extraction of pairing gap from the experimental level densities [l^ 
confirms that the pairing gap does not vanish at Tc but decreases as T increases, in line with 
the predictions by these approaches. 

The above mentioned approaches are based on the independent quasiparticles, whose 
occupation numbers follow the Fermi-Dirac distribution of free fermions. Dynamic effects 
such as those due to coupling to small-amplitude vibrations within the random-phase ap- 
proximation (RPA) are ignored. These effects have recently been explored by extending the 
self-consistent particle-particle RPA (SCRPA) to finite temperature using the double-time 



Green's function method 



17l |. Since the SCRPA fails in the region of strong pairing, where 



it should be replaced by the quasiparticle representation, it is highly desirable to develop a 
self-consistent quasiparticle RPA (SCQRPA) at finite temperature, which is workable with 
any value of pairing interaction parameter G. 



Recently, we have developed in Ref. 



18| a SCQRPA for the multilevel pairing Hamilto- 



nian and applied it to the Richardson model 19| at zero temperature. The derivation of the 
SCQRPA is based on a set of renormalized BCS equations, which include the corrections 
due to the QNF and the SCQRPA. The latter arise from the expectation values (AjA^j,) and 
{AjAj') in the correlated ground state. Here Aj is the product of two time-reversal con- 
jugated quasiparticle operators, a] and a^j, corresponding to the j-th orbital. Within the 
particle-particle (pp) SCRPA [20|], these expectation values overscreen the attractive pairing 
interaction, turning it into repulsion in agreement with the trend of the exact solutions of 
the Richardson model. For this reason, the expectation values (A^jA^j,) and {A^jAj') are 
called the screening factors. The goal of the present study is to extend the SCQRPA in Ref. 
jisl ] to non-zero temperature to explore the effects due to QNF as well as coupling to QRPA 
vibrations on the pairing properties of finite systems in a self-consistent way. 

The article is organized as follows. The derivation of the equations for quasiparticle 
propagation, which include the effects of QNF and SCQRPA corrections as well as cou- 
pling of quasiparticles to pair vibrations at finite temperature is presented Section [Tll Two 
approximation schemes will be considered, which are based on the thermal quasiparticle 
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representation without and including dynamic coupling to SCQRPA quasiparticle-pair vi- 
brations. In Section [Int the developed approach undergoes a thorough numerical test within 
the Richardson model as well as in realistic nuclei ^^Fe and ^^°Sn. The last section summa- 
rizes the article, where conclusions are drawn. 



II. FORMALISM 



A. Quasiparticle Hamiltonian 

The pairing Hamiltonian 

jm jj' mm'>0 

describes a set of particles with single-particle energies e^, which are generated by particle 
creation operators a^^ on j-th orbitals with shell degeneracies 2Qj (Qj = j + 1/2), and inter- 
acting via a monopole-pairing force with a constant parameter G. The symbol ~ denotes 
the time-reversal operator, namely = {—y~"^0'j-m- In general, for a two-component 
system with Z protons and neutrons, the sums in Eq. ([1]) run over all jriTir, j'r^'r, and 
Gr with r = {Z, N). This general notation is omitted here as the calculations in the present 
article are carried out only for one type of particles. 

By using the Bogoliubov's transformation from the particle operators, aj^ and ajm, to 
the quasiparticle ones, aj^ and ajm, 



the pairing Hamiltonian ([T]) is transformed into the quasiparticle Hamiltonian as follows 13|, 



3 j jj' jj' 

jj' jj' 

where Mj is the quasiparticle-number operator, whereas and Aj are the creation and 
destruction operators of a pair of time-reversal conjugated quasiparticles: 

^j = Yl ^\yn(^i^ = Y^^lrnOi^"" + a]_^aj_^) , (4) 

m=—Ui m=l 



A = [«] ® «]] = ^ E ' A- = ■ (5) 



' ''^ JU, ^ 

^ V J m=l 

They obey the following commutation relations 



[A, , A],] = 5,,,V, , where V, = l-^, (6) 



, A],] = 25,,. A], , [AT, , A,] = -25,,. Ay . (7) 

The functionals a, bj, Cj, djj', gj{j'), hjj', qjjf in Eq. ([3]) are given in terms of the coefficients 
Uj, Vj of the Bogoliubov's transformation, and the single particle energies ej as (See Eqs. 



(7) -(13) ofRef. [lag, e.g.) 

a = 2 5^ n,e,v^ -GiY, ^jU^v,)' - G Q.v^ , (8) 



bj = ^jiu] - v]) + 2GujVj Qj.Uj.Vj. + Gv^ , (9) 

i' 

Cj = 2^/fijejUjVj — Gy'Tlj^u'^ — f|) Qj/UjiVji — 2G^/fljUjVj , (10) 

j' 

^jj' = —G^/^lj^lj'{u'^Uy + VjVj,) = djtj , (11) 

9j{3) = GujVj^9~,{u], - v]) , (12) 

Q 

hjj' = -y^^j^j'iu'^jv], + v]u]) = hj,j , (13) 

Qjj' = -GujVjUjiVj, = qj>j . (14) 

By setting Qj = 1 in Eqs. ([H]) - dHD; one recovers the expressions for the case with Q 



doubly- folded levels of the Richardson model (See, e.g., Eqs. (12) - (18) of Ref. [18l]). 



B. Gap and number equations 

The derivation of the equation for the pairing gap that include the effect of correlations 



in the ground state has been presented briefly in Ref. [18(] for the Richardson model. For 
the clarity of the extension to finite temperature T, we give below the detailed derivation 
of the gap equation, which is applied to the more general quasiparticle Hamiltonian (IHl) and 
valid for T 7^ 0. 
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The coefficients uj and Vj of the Bogohubov's transformation ([2]) are determined by 
using the variational procedure, which minimizes the expectation value of the Hamiltonian 



Ti = H — XN in the grand canonical ensemble. This leads to the variational equations 21 1 



a(H)^cgO|H..(|„,.,).0, (15) 



duj dvj duj 
where (O) denotes the ensemble average of the operator O, 

The commutation relation [H, A'j] is found by using Eqs. and (171) as 

[n, A]] = 2b'^A] + |c; + J2 [djMf + 9Am' + h,rM 

+ 2 E { ^:'-(^") + + + A/;-' A^] I + E • (17) 

j' ^ ^ j' 

The ensemble average of the commutation relation ffTTj) is then given as 

where the functionals b'j and c'j are 

b'j = bj - A(m2 - v]) , c'j = Cj - 2\^jUjVj , (19) 

i.e. they have the same form as that of bj in Eq. (Q, and cj in Eq. (ITOl) . but with ej — A 
replacing ej at the right-hand sides. Inserting the explicit expressions for the functionals c'j 
from Eq. f|T9l) as well as gj{j') and 5'j'(j) from Eq. f|T2|) into the right-hand side of Eq. f|T8l) . 
and equalizing the obtained result to zero as required by the variational procedure (fT5l) . we 
come to the following equation, which is formally identical to the BCS one: 

2(e; - Gv] - \)u,vj - A,iu'^ - v]) = , (20) 

where, however, the single-particle energies e'^ are renormalized as 

G_ 

~ ^ V^A-D, 

The pairing gap is found as the solution of the following equation 



- ^AV(4 - ^1') {}^'A') + (A' A')) • (21) 



A ^ 



, (22) 
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which is level- dependent. The coefficients uj and Vj of the Bogoliubov's transformation ([2]) 
are derived in a standard way from Eq. fl20|) and the unitarity constraint u'j + Vj = 1. They 

where i?j are the quasiparticle energies 



E, = ^ie'^-Gv]-X)^ + ^ . (24) 

The particle-number equation is obtained by transforming the particle-number operator 
^ a}j^<ijm into the quasiparticle presentation using the Bogoliubov's transformation 
(|2]) and taking the ensemble average. The result is 



N = 2Y,^\v%V,) + \{l-{V,)) 



(25) 



The pairing gap Aj and chemical potential A, which is the Lagrangian multiplier in the 
variational equations (fT5|) . are determined as solutions of Eqs. ( 122!) and (!25|) . 

The right-hand side of Eq. fl22|) contains the expectation values {'D/Dji), whose exact 



treatment is not possible as it involves an infinite boson expansion series 

n 

article, following the treatment on Ref. p^, we use the exact relation 



22l | . In the present 



c \ r 

{Vp,,) = (V,) {Vy) + — ^ , with <5A/;-,, = (AT^.,) - {N,) {Ny) , (26) 

u L J u L J I 

and the mean-field contraction for the term bM^y 

5Mjf - 2njSf^Sjj' , SAff = rijil - rij) , (27) 
with the quasiparticle occupation number nj 

% = ^ = ^(1 - i^^)) , (28) 

to rewrite the gap equation fl22l) as a sum of a level-independent part, A, and a level- 
dependent part, SAj, namely 

Aj = A + 6Aj , (29) 

where 

5 J\f'^ 

A = G X] n,, {V,,)uyv,, , 6 A, = ^G^^ujv, . (30) 



The quantity SMf in Eqs. ( I27I) and (|30D . is nothing but the standard expression for the QNF 



corresponding to the j-th orbital 



^. Using Eqs. fl23l) and fl30|) . after simple algebras, 



we rewrite the gap (1291) in the following form 

A,- = ^ V n,,iV,,)^ , where G,- = g(i~ G-^^^ 



(31) 



C. Finite-temperature BCS with quasiparticle number fluctuations 

1. Without particle-number projection (FTBCSl) 

The gap equation fl3Tl) is remarkable as it shows that the QNF renormalizes the pair- 
ing interaction G to Gj. The conventional finite-temperature BCS (FTBCS) gap equation 
Aj = A is recovered from Eq. ( 13T1) when the following assumptions simultaneously hold: 

i) Independent quasiparticles: rij = nj^, where is the Fermi-Dirac distribution of 
non-interacting fermions 

(32) 



n 



FD 



1 



^1 ^(3E, + I ' 

ii) No quasiparticle number fluctuation: 6Afj = , 

iii) No screening factors: (^^UJ,) = {A]Af) = in Eq. (|2T1). 

These three assumptions guaranty a thermal quasiparticle mean field, in which quasiparti- 
cles are moving independently without any perturbation caused by the QNF and/or coupling 
to multiple quasiparticle configurations beyond the quasiparticle mean field. Among these 
configurations, the simplest ones are the small-amplitude vibrations (QRPA corrections). 
From these assumptions, one can infer that releasing assumption ii) allows us to include the 
effect of QNF, provided the quantal effect of coupling to QRPA vibrations is negligible, i.e. 
assumption iii) still holds. In the present article, this approximation scheme, for which i) 
and iii) hold, whereas 6Afj ^ 0, is referred to as the FTBCSl. 



1 The definition ^ for the QNF bM] is different from that in Eq. (32) of Ref. [jj] by a factor 2 as this 
factor is now put in front of to have the complete shell degeneracy 20. j . 
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2. With Lipkin-Nogami particle-number projection (FTLNl) 



The problem of particle-number violation within the BCS theory is usually resolved in 
the simplest way by means of an approximated particle-number projection (PNP) before 
variation called the Lipkin-Nogami (LN) method 23|]. In Ref. [18| this method has been 
applied to the BCSl and the resulting approach is called the LNl. For the case with Vtj ^ 
1 and level-dependent gap fl29|) at T 7^ 0, the corresponding finite-temperature LNl 
equations have the form 



A, 



N 



2 ^jPj 



e'j + {AX2-G)v^ , X = Xi + 2X2{N + 1) , (33) 



where 



2 SAT? ^ ^ 




The coefficient A2 is given as 



18| 




311 ' 



(e, - xy + K] . 



(34) 
(35) 



I Til 



(36) 



where Tj = Tjj. This FTBCSl including the approximated PNP within the LN method 
is referred to as FTLNl in the present article. It is worth pointing out that, being an 
approximated projection that corrects for the quantal fluctuations of particle number within 
the BCS theory, the LN method in the present formulation is not sufficient to account for 
the thermal fluctuations (QNF) around the phase transition point T ~ T^. as well as at high 
T . Another well-known defect of the LN method is that it produces a large pairing gap 
(pairing correlation energy) even in closed-shell nuclei, where there should be no pairing 
gap. The source of this pathological behavior is assigned to the fast change of A2 at the 
shell closure, which invalidates the truncation of the expansion at second order 241. In Ref. 

n 

|15l |. it has been demonstrated within the MBCS theory that the projection-after- variation 
(PAV) method offers much better results, which are closer to the exact solutions. The PAV 
at T 7^ 0, however, is much more complicated than the LN method. Therefore, we prefer to 
devote a separate study to its application to the BCSl. 
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D. Finite-temperature BCS with quasiparticle-number fluctuation and dynamic 
coupling to SCQRPA vibrations (FTBCSl+SCQRPA and FTLNl+SCQRPA) 



As has been mentioned in the preceding section, within the quasiparticle mean field, 
the expectation values {AjA^j,) and {A^jAj') at the right-hand side of Eq. ([5T]) are always 
zero [Assumption iii)]. They cannot be factorized into the products of expectation values 
of quasiparticle-number operators within the thermal quasiparticle mean field because such 
crude contraction is tantamount to artificially breaking the pair correlators ([5]) (See the 
Appendix A). Therefore, to account for the correlations beyond the quasiparticle mean 
field, these screening factors should be estimated, at least, within the SCQRPA, where they 
can be expressed below in terms of the forward- and backward going amplitudes, X^" and 
y!^, of the SCQRPA operators (phonons) as 

^I^J:{^a]-^a,), ,37, 



The renormalization factor a/ {'Dj) is introduced in Eq. fl37|) to ensure that the SCQRPA 
operators Q^^ and remain bosons within the thermal average (fT6|) . preserving the exact 
commutation relation ([6]). This leads to the orthogonality relation for the Xj' and yj 
amplitudes in the conventional form as 

Y: {x^xf - ypf) = s,,' , (38) 

j 

which can be easily verified by calculating {[Q^, Qj^/]) and requiring that the result to be 
equal to 6^^/. The inverse transformation of Eq. fl37|l reads 



provided the following conventional closure relations hold 



In general, operator Qj^{JM) at T 7^ also contains the terms ^ B^^, (JM) = [aj (g) aj']^^ and Bjj' [JM) 
apart from those with A]j,{JM) and Ajj'UM) because of the relation {[Bjji{JM),B\^.,^{J' M')]) = 
5j.j'5MM'5jj^5jiji^ {rij — Tij') ^ for j ^ j' [25|, |26|. In the present article, where J = M = 0, and hence 
j = j', this relation vanishes. 
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1. Screening factors 



Using the inverse transformation (139|) . we obtain the expectation values {J^-J^-,) and 



r J' 



{A] Ay) at T 7^ in the form 

{A\A,.) 



{A]A]) 



Y.y'^r + E ( U^r^Ql') + ^jI^'(2iQ.') j , (42) 

(43) 



where the following shorthand notations are used 



^jr = ^'^r + ypf ' = -^/^j' + yf^r ■ 

taking into account the symmetry property (Qj^Q^^,) = {Q^Qf^'). Using now the definition 
( !37ll . we express the expectation values {Qjj^Q^') and {Qj^Q^^,) in terms of Xjj> (i.e. {A^jAj')), 
Ujji (i.e. (^j^L)), and amplitudes and as 



where 



(44) 
(45) 

(46) 



Inserting Eqs. ( l44l) and ( 1451) into the right-hand sides of Eqs. ( 1411) and (1421) . after some 
simple algebras, we obtain the following set of exact equations for the screening factors (14T|) 
and (in 



= E^^ + E - ^''^o > (47) 



= E3T'^+E3^'(^'3^'-f^''^')- 



(48) 
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The derivation of the SCQRPA equations at finite temperature is proceeded in the same 
way as has been done at T = 0, and is formally identical to Eqs. (46), (56), and (57) of 
Ref. jl^ so we do not repeat them here. Notice that the expectation values {VjVji) in the 
submatrices A and B in Eqs. (56) and (57) of Ref. [isl] are now calculated by using Eqs. fl2^ 
and (127|1 . The approach that solves the number and gap equations ( 1251) . (129]) - (1271) . as well 
as equations for the screening factors (H71) and (HHl) selfconsistently with the SCQRPA ones 
at T 7^ 0, where all the assumptions i) - iii) cease to hold, is called the FTBCSl+SCQRPA 
in the present article. The corresponding approach that includes also PNP within the LN 
method is called as FTLNl+SCQRPA. 



2. Quasiparticle occupation number 

To complete the set of FTBCSl+SCQRPA equations we still need an equation for the 
quasiparticle occupation number nj defined in Eq. ( l28l) . Here comes the principal difference 
of the FTBCSl+SCQRPA compared to the zero-temperature SCQRPA since nj should be 
calculated selfconsistently from the SCQRPA taking into account dynamic coupling between 
quasiparticles and SCQRPA phonons at T 7^ in an infinite hierarchy of algebraic equations. 
The quasiparticle propagator found as the formal solution of this hierarchy of equations is 
different from that for free quasiparticles by the mass operator, which reflects the effects of 
coupling to complex conflgurations. Since the latter cannot be treated exactly, approxima- 



tions have to be made to close the hierarchy. Following the same line as in Ref. 17|, we 



e oc- 



as]. 



derive in this section a set of equations for the quasiparticle propagator and quasipartic 
cupation number rij at T 7^ by using the method of double-time Green's functions [27 
To close the hierarchy of equations, we lower the order of double-time Green's functions by 
ying the standard decoupling approximation introduced by Bogoliubov and Tyablikov 



app] 


yin 


27I, 


28] 



By noticing that the only term in the quasiparticle Hamiltonian ([3]) that cannot be taken 
into account within either the BCS theory or the SCQRPA is the sum containing gj{j') 
functionals, we effectively rewrite Hamiltonian Ti in Eq. (|T5|) as 

^^// = E(^^- + E In'^r)^^ + E ^/^^i^M + E Vf^fM + 2m) • (49) 
j 3' M J> 

The flrst sum at the right-hand side of this representation describes the part of the quasipar- 
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tide Hamiltonian ([3]), which cannot be expressed in terms of phonon operators (1371) . Within 
the BCS theory, where the part containing qjji does not contribute whereas the term ~ Gv'j 
and the QNF are neglected, one obtains h'^ = Ej. In this case, this sum corresponds to the 
quasiparticle mean field. The second sum describes the SCQRPA Hamiltonian after solving 
the SCQRPA equations, which give the amplitudes Xj", y^, and the SCQRPA energies 



UJ 



The last sum represents the coupling between the quasiparticle and phonon fields, which is 
left out from the BCS (FTBCSl) and the QRPA (SCQRPA). This sum is rewritten here in 
terms of Afj and SCQRPA operators by using the inverse transformation (139|) . The vertex 



Vj^ obtained after this transformation has the form 



vr = Y.9Af)^{i^f)i^i^+y^')- (50) 
j' 

Given that Afj commutes with Qj^ within the SCQRPA, such effective representation of 
the quasiparticle Hamiltonian causes no double counting between the first two sums at the 
right-hand side of Eq. (1491) . but becomes convenient for the derivation of the quasiparticle 
Green's function, which includes the coupling to SCQRPA modes, because the first sum 
is activated only in the quasiparticle space, whereas the second sum functions only in the 
phonon space. 



28|, we introduce the 



Following closely the procedure described in Section 8.1 of Ref. 
double-time retarded Green's functions, which describe 

a) The quasiparticle propagation: 

G,{t-t') = {{a,{ty,al{t'))) , (51) 

b) Quasiparticle-phonon coupling: 

T.^{t-t') = {{a,{t)Q,{ty,a]{t'))) , r+ (t - t') = {{a,{t)Ql{t); a]{t'))) . (52) 

The magnetic quantum number m in aj^ and ajm is omitted hereafter for simplicity as the 
results below do not depend on m. The definitions (15T1) and fl52l) use the standard notation 
G^{t - t') = {{A{t); B{t'))) = -i9{t - t'){[A{t), B{t')]) for the double-time retarded Green's 
function Gr{t — t') built from operators A{t) at time t and B(t') at time t' . The advantage 
of using the double-time retarded Green's function is that this type of Green's function can 
be analytically continued into the complex energy plane. The imaginary part of the mass 
operator in this analytic continuation corresponds to the quasiparticle damping caused by 
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the quasiparticle-phonon coupling. This method is free from any constraints of perturbation 
theory. 

Applying the standard method of deriving the equation of motion for the double-time 
Green's function, namely 
3C (f — f') 

^ ^ = 6{t-t'){[A{t),B{t)]^) + {{[A{t),H{t)]-B{t'))) , (53) 

to the Green's functions ( ISTl) and ( l52|) with the effective Hamiltonian (l49ll . we find for them 
a set of three exact equations 

^ %^ ^ = Sit - t') + E,G,{t - t') + E [^U^ - + - ^')] , (54) 

dT~ it — t') ~ 

^^^^ = {E, + u:,)rj^{t - t') + E Vf'(Mt) Kit) + Q,,{t)] Q,{t)- a]{t'))) 

+ Y.^,^^{a,{t)U,,{t)-a\t'))) , (55) 

r 

QY+ it - t') ~ 

^^^^ = iE,-uj,)rl^{t-t') + Y,Vf{{a,{t)[Ql{t) + QAt)]Qiity,a]im 

-J2^p{MtWAty,c^\t'))), (56) 

i' 

where 

E, = 6; + g,, . (57) 

The last two equations, Eqs. (l55l) and (l56l) . from this set contain higher-order Green's 
functions, which should be decoupled so that the set can be closed. Following the method 
proposed by Bogoliubov and Tyablikov 27|, we decouple the higher-order Green's functions 
at the right-hand side of Eqs. flS^ and fl56l) by pairing off operators referring to the same 
time, namely 

(Mt) [QlA) + QAt)] QM <^]it'))) ^ s,,,u,G,it - A , 

(Mt) Kit) + QAt)] Qiit); a]{A)) ^ 6,A^ + '^,)GA " t') , (58) 
{{a.mAty, a\A)) ^ (1 - n,)GA - A ■ 
As the result of this decoupling, Eqs. (!55|) and (l56ll become 
QYT (t - 1') 

^^^^ = iEj + oo,)r-A - A + yfi^ - n, + v,)GA - A , (59) 
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{E, - uj,)TUt - t') + Vf{n, + u,)G,{t - t') 



(60) 



Taking the the Fourier transforms of Eqs. ( l54l) . ( l59l) . and ( l60l) into the (complex) energy 
variable E, one obtains three equations for three Green's functions Gj{E), TJ^{E), and 
Tj'^(E). Eliminating two functions T^^{E) by expressing them in terms of Gj{E) and in- 
serting the results obtained into the equation for Gj{E), we find the final equation for the 
quasiparticle Green's function Gj{E) in the form 

1 1 



G,{E) 



"27^ E - Ej - Mj{E) 



where the mass operator Mj{E) is given as 



E — Ej — uj^ 



+ 



E-Ej+uj^,, 



(61) 



(62) 



In the complex energy plane E = uo ±ie {uj real), the mass operator (162|) can be written as 



where 



Mj{uj ± ie) = Mj{u) =F ilj{oo) , 



;i - Hj + v^){u - Ej - u^) ^ {rij + v^){u - Ej + u^) 



{uo - Ej - uo^y + £2 



{uo - E, + u^f + 52 



7.M = ^E(^') 



l-rij^ Uf, 



rij + z/^ 



{UJ - E, - UJ^f + £2 (c^ _ + cu^)2 + ^2 



The spectral intensity Jjiuj) of quasiparticles is found from the relation 

Gj{u + ie) - Gj^cj - ie) = -iJj{uj){e^'' + 1) , 

and has the final form as 27|, |28 1 

J M = 1 7,(^)(e^" + l)-^ 

vr[^_i5^._M,(c.)P + 72(^) • 

Using Eq. (jSZ]), we find the quasiparticle occupation number rij as the limit t = t' of the 
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(64) 
(65) 

(66) 

(67) 



correlation function 



The final result reads 
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(68) 



(69) 



In the limit of small quasiparticle damping 7^(0;) 0, the spectral intensity Jj{uj) becomes 
a 5-function, and rij can be approximated with the Fermi-Dirac distribution [exp(/3£^j) + l]^^ 
at 'jjiE'j) , where Ej is the solution of the equation for the pole of the quasiparticle 
Green's function Gj{uj), namely E'j — Ej — Mj{E'-) = , whereas the quasiparticle damping 
at = Ej due to quasiparticle-phonon coupling is given by •jj^Ej). 

We have derived a closed set of Eqs. (IMI) . (l65il . and (169|) for the energy shift Mj, damping 
7j, and occupation number rij of quasiparticles. which should be solved self-consistently with 
the SCQRPA equations at T 7^ with the screening factors calculated from Eqs. (1471) and 
( l48l) . The quasiparticle occupation number rij obtained in this way is used to determine the 
pairing gap from Eq. fl29l) . These equations form the complete set of the FTBCSl+SCQRPA 
equations for the pairing Hamiltonian ([T]), where the dynamic effect of quasiparticle-phonon 
coupling is self-consistently taken into account in the calculation of quasiparticle occupation 
numbers. 



III. ANALYSIS OF NUMERICAL RESULTS 
A. Ingredients of calculations 

We test the developed approach by carrying out numerical calculations within a schematic 
model as well as realistic single-particle spectra. For the schematic model, we employ the 
Richardson model having Q doubly-folded equidistant levels with the number Q of levels 
equal to that of particles, N. This particle-hole symmetric case is called the half-filled one 
as in the absence of the pairing interaction {G = 0), all the lowest Q/2 levels are occupied 
by N particles with 2 particles on each level. The level distance is taken to be 1 MeV to 
have the single particle energies ej = j MeV with j = 1, . . . ,Q. The results of calculations 
carried out within the FTBCSl, FTLNl, FTBCSl+SCQRPA, and FTLNl+SCQRPA at 



various and G will be analyzed. As this model can be solved exactly [29|, for the sake 
of an illustrative example, we will compare the predictions by these approximations with 
the exact results obtained for = 10 and G = 0.4 MeV after extending the latter to finite 
temperature. Such extension is carried out by averaging the exact eigenvalues over the 



canonical ensemble of A^ particles 14 1. 



For the test in realistic nuclei, Fe and Sn, the neutron single-particle spectra for the 
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bound states are obtained within the Woods-Saxon potentials at T = 0, and kept unchanged 
as T varies. The parameters of the Woods-Saxon potential for ^^°Sn take the following values: 
V = -42.5 MeV, V^i^ = 16.7 MeV, a = = 0.7 fm, R = 6.64 fm, and Ri, = 6.46 fm. The 
full neutron spectrum for ^^°Sn spans an energy interval from around —37 to 7.5 MeV for 
^^°Sn. From this spectrum the calculations use all 22 bound orbitals with the top bound 
orbital, 1^13/2, at energy of —0.478. For ^^Fe, as we would like to compare the results of our 
approach with the predictions by the finite-temperature quantum Monte Carlo (FTQMC) 



SOj , the same single-particle energies from Table 1 of Ref . 



3011 for 



method reported in Ref. 
^^Fe and the same values for G therein are used in calculations. Given the large number of 



results reported in Ref. [30||, we choose to show here only one illustrative example for the 
pf shell. 



The main quantities under study in the numerical analysis are the level-weighted gap 

A = ^^^-^ , (70) 



total energy S = (H), and heat capacity C = dS/dT. By using PNP within the LN method, 
the internal energy has an additional term due to particle- number fluctuations AA^^ {23!, 
namely 

Within the FTLNl, the particle-number fluctuations AA^^ consist of the quantal fluctuation, 
AiVqp, and statistical one, AA^|p, which are calculated following Eqs. (16) and (17) in Ref. 



3l| . respectively. Within the FTLNH-SCQRPA, a term 6Nsc due to the screening factors 



should be added, so that 

AN^ = AN^^+ANIp+6Nsc , SNsc = 8 ^ ,/Q^ujVjUj>vy [{A]A],) + {A]Af)] . (72) 

jj' 

The integration in Eq. ( l69il is carried out within the energy interval — cul < < ci^l with 
cul = 100 MeV and a mesh point Auj < 0.02 MeV. Since the integration limit is finite, 
the integral flB^ is normalized by J^^^ Jj{uj)[exp{Pu) + l]du. The results obtained within 
the FTBCSl+SCQRPA (FTLNl+SCQRPA) by using a smearing parameter e < 0.2 MeV 
[in calculating the mass operator (!64|) and quasiparticle damping (!65|) ] are analyzed. They 
remain practically the same with varying e up to around 0.5 MeV. 
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T (MeV) T (MeV) 

FIG. 1: (Color on line) Level-weighted pairing gaps A obtained within the FTBCSl as functions 
of temperature T at various values of pairing parameter G (in MeV) indicated by the figures near 
the lines for several values of particle number N . Open circles on the axes of abscissas in panels 
(a) and (b) mark the values Ti of temperature, where the FTBCSl gap turns finite at low G. Full 
circles denote temperature Tc, where the gap vanishes, and T2, where it reappears. 

B. Results within Richardson model 

1. Effect of quasiparticle-number fluctuation 

It is well known that, below a critical value Gc of the pairing interaction parameter, the 
conventional BCS theory has only a trivial solution (A = 0). At G > Gc, the FTBCS gap 
decreases with increasing T up to a critical value of T = T^, where it collapses, and the 
system undergoes a sharp SN-phase transition. The behavior of the pairing gap within the 
FTBCSl theory can be inferred from Eq. fl3Tl) . As a matter of fact, the increase of the 
QNF 6Af^ with T leads to an increase of Gj, whose consequences are qualitatively different 
depending on the magnitude of G and particle number N. These features can be seen in Fig. 
[H where the level- weighted pairing gaps A obtained within the FTBCSl theory at various 
values of the pairing interaction parameter G for several particle numbers are displayed as 
functions of temperature T. They can be classified in three regions below. 
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In the region of strong coupling, G ^ Gc, where the BCS equations have non-trivial 
solutions at T = 0, and is sufficiently large so that Gj ^ G, the gap Aj in Eq. fl3T|) 
never collapses since whenever T reaches the value T^. where the BCS gap obtained with 
parameter G collapses, the gap is always positive given Gj ^ G with a renormalized 
critical temperature ^ T^. In this way, the sharp SN-phase transition never occurs as Aj 
remains always finite at Tc < T ^ Tc with Tc continuously becoming larger with T. If G is 
sufficiently large the QNF may become so large at high T that the level-dependent part 6Aj 
in Eqs. (1291) and (!30l) starts to dominate and the total gap A^ will even increase with T. 
This effect is stronger when the particle number is smaller. As seen in Fig. [1], in contrast to 
the FTBCS gap, which collapses at T^, the FTBCSl gaps shown as the thick solid lines are 
always finite. For > 6, the gaps decrease monotonously as T increases up to T = 4 MeV. 
This feature qualitatively agrees with the findings within alternative approaches to thermal 
fiuctuations mentioned in the Introduction. 

In the region of weak coupling, G < Gc, where the pairing gap is zero at T = 0, the 
increase of Gj with T makes it becomes significantly greater than Gc at a certain T = Ti, 
allowing a non-trivial solution of the gap equation. This feature is demonstrated by the 
dotted lines in Figs. [1] (a) and [1] (b), where Ti (> 2 MeV) is marked by an open circle. 
Since the difference between the FTBCSl gap Aj and the conventional FTBCS one. A, is 
the gap 6Aj in Eqs. fl29|) and fl30|) . which arises because of the QNF it is obvious that 
the finite gap at T > Ti is assisted by the QNF. 

In the transitional region, where G is slightly larger than Gc, it may happens that, 
although increases with T, it is still too small so that Gj is only slightly larger than 
G, and so is Tc compared to Tc. As a result, the gap collapses at T = Tc which is slightly 
larger than Tc. As T increases further, the mechanism of the weak-coupling region is in 
effect, which leads to the reappearance of the gap at T = T2 > Tc. In Fig. [1], these values 
Tc and T2 are denoted by full circles on the axes of absiccas for the cases with A^ = 6, 10, 
20, 50 with G = 0.6, 0.4, 0.3, and 0.24 MeV, respectively. With increasing G, it is seen that 
Tc increases whereas T2 decreases so that at a certain G these two temperatures coalesce. 
The value Gm where f c = T2 = Tm is found to be 0.651, 0.51815, 0.40205 and 0.3095 MeV 
for A^ = 6, 10, 20 and 50, respectively, i.e. decreases with increasing A^. The gap obtained 
with G = Gm is seen decreasing with increasing T from to Tm, where it becomes zero. 
Starting from Tm the gap increases again with T. The value Tm is found increasing with T 
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FIG. 2: (Color on line) Level-dependent pairing gap Aj ([29|) and level-weighted pairing gap A (fTOl) 
obtained within the FTBCSl as functions of temperature T for = 20 and G = 0.44 MeV. Thick 
solid lines represent the level-weighted gaps A. Thin solid lines denote the level-dependent gaps 
Aj corresponding to the j-th orbitals, whose level numbers j are marked at the lines. Dashed and 
dotted lines stand for the level-independent part (quantal component), A, and the level-dependent 
one (thermal component), SAj, of the FTBCSl gap Aj respectively. 

from Tm ~ 1.2 MeV for = 6 to Tm ~ 1.7 MeV for = 50. At G > Gm the gap remains 
finite at any value of T. For small A^, the strong QNF even leads to an increase of the gap 
with T at high T as seen in the cases with A^ = 6, and Gm < G < 1.2 MeV. With increasing 
A^ the high-T tail of the gap gets depleted, showing how the QNF weakens at large A^. 

The curious behavior of the level-weighted gap at weak coupling, where it appears at a 
certain T = Ti, and in the transitional region, where it collapses at and reappears at 
T2, may have been caused by the well-known inadequacy of the BCS approximation (and 



BCS-based approaches) for weak pairing [32]. Even at T = 0, Ref. [29| has shown that, 
whereas the exact solution predicts a condensation energy of almost 2 MeV in the doubly- 
closed shell ^'^Ca, the BCS gives a normal Fermi-gas solution with zero pairing energy. It is 



expected that a proper PNP such as the number-projected HFB approach in Ref. |33j, if it 
can be practically extended to T 7^ 0, will eventually smooth out the transition points Ti as 
well as Tc and T2 in Fig. [H 

To have an insight into the source that causes the high-T tail of the FTBCSl gap we plot 
in Fig. [2] the examples for the level-weighted gaps A fITOl) along with the level-dependent 



gaps Aj ( !29l) . which are obtained for A^ = 20 and G = 0.44 MeV. It is seen from this figure 
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that the level-independent part (quantal component) A of the gap [dashed lines in Fig. 
[2] (b)] also has a high-T tail although it is much depleted compared to the total gap Aj, 
which includes the level-dependent part SAj. This figure also reveals that the QNF has the 
strongest effect on the levels closest to the Fermi surface, which are the 10th and 11th levels. 
In this figure, the results for the 11th level are not showed as they coincide with those for 
the 10th one due to the particle-hole symmetry, which is well preserved within the FTBCSl. 
For the rest of levels, the effect of QNF is much weaker. With increasing the particle number 
A^, the number of levels away from the Fermi surface becomes larger, whose contribution in 
the gap A outweighs that of the levels closest to the Fermi surface. This explains why the 
high-T tail of the level-weighted gap A is depleted at large A^. When becomes very large, 
this tail practically vanishes as the total effect of QNF becomes negligible. In this limit, the 
temperature dependence of the pairing gap approaches that predicted by the standard BCS 
theory, which is well valid for infinite systems. 

2. Corrections due to particle-number projection and SCQRPA 

Show in Fig. [3] are the level- weighted pairing gaps A, total energies S, and heat 
capacities C, obtained within the FTBCS, FTBCSl, FTLNl, FTBCSl+SCQRPA, and 
FTLNl+SCQRPA for the systems with A^ = 10 (G = 0.4 MeV) and A^ = 50 (G = 0.3 
MeV). As we want to see the effect of QNF for the case with small A(T = 0) without 
any phase transition points at and T2, we choose to neglect, for this particular test, the 
self-energy term — GfJ in the single-particle energy. For A^ = 10 e.g., this increases the 
gap at T = by around 14%, to around 0.8 MeV, but the change in the total energy is 
found to be negligible. Different from the common practice, which usually neglects the terms 
~ —G J2j ^j'^j ill calculating the total energy S, the latter is calculated in the present article 
by averaging the complete pairing Hamiltonian ([3]). For A^ = 10 and G = 0.4 MeV e.g., this 
causes a shift of total energy down by around 2 MeV (~ 8%) and 1 MeV (~ 10.4%) at T = 
and 4 MeV, respectively. 

As has been discussed in Sec. Ill C 2\ Fig. [3] demonstrates that, although the LN method 
significantly improves the agreement between the predictions by the FTBCSl theory with 
the exact results for the pairing gap and total energy at low T, it fails to do so at T > Tc, 
where all approximated results for the pairing gap coalesce and clearly differ from the exact 
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FIG. 3: (Color on line) Level- weighted pairing gaps (a, d), total energies (b, e), and heat capacities 
(c, f) as functions of temperature T, obtained for = 10 [(a) - (c)], and = 50 [(d) - (f)]. The 
dotted, thin sohd, thick soUd hues show the FTBCS, FTBCSl, and FTBCSl+SCQRPA results, 
respectively. The predictions by the FTLNl and FTLNl+SCQRPA are presented by the thin and 
thick dashed lines, respectively. The dash-dotted lines in (a) - (c) denote the exact results. The 
calculations of the mass operator and quasiparticle damping within the SCQRPA were performed 
using £ = 0.05 MeV. 

result (for = 10). The reason is partly due to the fact that, strictly speaking, there is no 



pairing gap in the exact solution [ij]. The dash-dotted line, representing the exact result in 
Fig. [3] (a) is the effective gap (canonical gap) extracted from the pairing energy. The latter 
is the difference between the exact total energy and the that of the single-particle mean field 
(Hartre-Fock) energies. The canonical gap includes correlations caused by the fluctuations 
of the order parameter, only a part of which is taken into account within the FTBCSl in 
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terms of QNF. It reduces to the BCS pairing gap only within the mean field approximation 
and the grand canonical ensemble 

The corrections caused by the SCQRPA are found to be significant for small (A^ = 10), 
in particular for the pairing gap in the region < T < 1.5 MeV [Fig. [3] (a)]. At T < T^, 
the predictions by the FTLNl+SCQRPA are closer to the exact results than those by the 
FTBCSl+SCQRPA. At T > Tc both approximations offer nearly the same results. They 
produce the total energies and heat capacities, which are much closer to the exact values 
as compared to the FTBCSl and FTLNl results, as shown in Figs. [3] (b) and [3] (c). What 
remarkable here is that the SCQRPA correction indeed smears out all the trace of the SN 
phase transition in the pairing gap as well as energy and heat capacity. 

For large A^ (A^ = 50), the effect of SCQRPA corrections is much smaller, although still 
visible. It depletes the spike, which is the signature of the SN phase transition around Tc 
in the heat capacity, leaving only a broad bump between < T < 2 MeV [Fig. [3] (f)]. 
The exact results are not available because, for large particle numbers, one faces technical 
problems of diagonalizing matrices of huge dimension, all the eigenvalues of which should be 
included in the partition function to describe correctly the total energy and heat capacity. 



C. Results by using realistic single-particle spectra 

The level-weighted gaps, total energies, and heat capacities, obtained for neutrons in ^^Fe 
and ^^°Sn within the same approximations are displayed as functions of T in Fig. |H The 
results of calculations for 10 neutrons in the I/7/2 2^3/22^1/21/5/2 shell using G = 25/26 MeV 
are plotted in Figs. H] (a) - H] (c) as functions of T within the same temperature interval 



as that in Ref. [3(|. They clearly show that the SCQRPA corrections bring the FTBCSl 
(FTLNl) +SCQRPA results closer to the predictions by the FTQMC method for the total 
energy and heat capacity (No results for the pairing gap are available within the FTQMC 



method in Ref. 30||). In heavy nuclei, such as ^^°Sn, the effects caused by the SCQRPA 
corrections are rather small on the pairing gap and total energy. In both nuclei, the pairing 
gaps do not collapse at T = T^, but monotonously decrease with increasing T, and the 
signature of the sharp SN-phase transition seen as a spike at T = Tc in the heat capacities 
is strongly smoothed out within the FTBCSl+SCQRPA. 
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FIG. 4: (Color on line) Level-weighted pairing gaps, total energies, and heat capacities for 10 
neutrons in the ^f-j/2^Pi/2'^Pi/2^fb/2 shell of ^^Fe and all neutron bound states of ^^'^Sn as functions 
of T (e = 0.1 MeV). Notations are as in Fig. [3l In (b) and (c), the predictions by the finite- 
temperature quantum Monte Carlo method [30] are shown as boxes and crosses with error bars 
connected by dash-dotted lines. 

D. Self-consistent and statistical treatments of quasiparticle occupation numbers 

The quasiparticle occupation numbers rij as predicted by the FTBCSl and FT- 
BCSl+SCQRPA for all quasiparticle levels in the system with = 10, G = 0.4 MeV, 
and for the orbitals within the (50 - 82) shell in i20Sn {G = 0.137 MeV) are shown in Fig. [5] 
as functions of T. While the ph symmetry is preserved within the FTBCSl {rij = n™) in the 
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FIG. 5: (Color on line) Quasiparticle occupation numbers for = 10 with G = 0.4 MeV (a, b) 
and ^^'^Sn with G = 0.137 MeV (c) as functions of T. In (a) and (b) the solid lines are predictions 
within FTBCSl+SCQRPA for the levels numerated by the numbers in the circles starting from the 
lowest ones. The dashed lines, numerated by the italic numbers, show the corresponding results 
obtained within the FTBCSl. In (c) predictions for the neutron orbitals of the (50 - 82) shell in 
^20Sn, obtained within the FTBCSl and FTBCSl+SCQRPA, are shown as the dashed and solid 
lines, respectively. 

sense that the values for raj^ are identical for the single-particle levels located symmetrically 
from the Fermi level [Compare the dashed lines in Figs. [5] (a) and [5] (b)], it is no longer the 
case after taking into account dynamic coupling to SCQRPA vibrations. This is particularly 
clear in light systems [See the solid lines in Figs. [5] (a) and [5] (b)]. This deviation of nj from 
the Fermi-Dirac distribution of free quasiparticles, however, turns out to be quite small in 
realistic heavy nuclei, such as ^^°Sn, as shown in Fig. [5] (c). 



E. Comparison between FTBCSl and MBCS 

In Refs. Q, Q, Q, Q the MBCS theory has been developed, which also produces 
a nonvanishing pairing gap at high T. Therefore, it is worthwhile to draw a comparison 
between the MBCS theory and the present one. Both approaches include the same QNF 
fl27|) as the microscopic source, which smoothes out the sharp SN-phase transition and leads 
to the high-T tail of the pairing gap. This high-T tail has been shown to be sensitive to the 
size of the configuration space in either approach. However, due to different assumptions in 
these two approaches, the functional dependences of 6Aj on the QNF 6J\fj are different. As 
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a result, the FTBCSl gap is level- dependent, whereas the MBCS one is not. The most im- 
portant advantage of the FTBCSl over the MBCS theory is that the solution of the FTBCSl 
gap equation (1291) is never negative. Moreover, at moderate and strong couplings, where 
the FTBCSl gap is finite, its behavior as a function of temperature bears no singularities 
in any configuration spaces for any value of > 2. The MBCS gap, on the other hand, 
is free from singularities only up to a certain temperature Tm, which is around 1.75 - 2.3 
MeV within the Richardson model with Q = N = 10 and increases almost linearly with N 
to reach Tm ~ 24 MeV hi n = N = 100 [14] (For detail discussions see Refs. Q, Q and 
references therein). However, the mean-field contraction used to factorize the QNF within 
the FTBCSl to the form flTTI) may have left out some higher-order fluctuations, which can 
enhance the total effect of the QNF. It might also be the reason that causes the phase tran- 
sition temperatures Ti, Tc and T2 at weak coupling and in the transitional region, discussed 
in Sec. IIII B 1[ Meanwhile, the MBCS theory is based on the strict requirement of restoring 
the unitarity relation for the generalized single-particle density matrix HI], which brings in 
the QNF 6Af^ (ETj) without the need of using a mean-field contraction. As a result, the effect 
of QNF within the MBCS theory is stronger than that predicted within the FTBCSl and/or 
FTBCSl+SCQRPA, which can be clearly seen by comparing, e.g.. Fig. |4](d) above and Fig. 



4 of Ref. 



111 ]. Whether this means that the secondary Bogoliubov's transformation prop- 



erly includes or exaggerates the effect of coupling to configurations beyond the quasiparticle 
mean field within the MBCS theory remains to be investigated. Another question is also 
open on whether the MBCS theory can be improved by coupling the modified quasiparticles 
to the modified QRPA vibrations. The answer to these issues may be a subject for future 
study. 



IV. CONCLUSIONS 



The present work extends the BCSl+SCQRPA theory, derived in Ref. [18| for a multilevel 
pairing model, to finite temperature. The resulting FTBCSl+SCQRPA theory includes the 
effect of QNF as well as dynamic coupling of quasiparticles to pairing vibrations. This theory 
also incorporates the corrections caused by the particle-number projection within the LN 
method. 

We have carried out a thorough test of the developed approach within the Richardson 
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model as well as two realistic nuclei, ^^Fe and ^^°Sn. The analysis of the obtained pairing 
gaps, total energies, and heat capacities leads to the following conclusions: 

1) The FTBCSl (with or without SCQRPA corrections) microscopically confirms that, in 
the region of moderate and strong couplings, the quasiparticle-number fluctuation smoothes 
out the sharp SN phase transition, predicted by the FTBCS theory. As a result, the gap 
does not collapse at T = Tc, but has a tail, which extends to high temperature T. 

2) The correction due to the particle-number projection within the LN method to the 
pairing gap is significant at T ^ Tc, which leads to a steeper temperature dependence of 
the pairing gap in the region around Tc. At the same time, the SCQRPA correction smears 
out the signature of a sharp SN phase transition even in heavy realistic nuclei such as ^^°Sn. 

3) The dynamic coupling to SCQRPA vibrations causes the deviation of the quasiparticle 
occupation number from the Fermi-Dirac distribution for non-interacting fermions. However, 
for a reahstic heavy nucleus such as ^^°Sn, this deviation is negligible. Consequently, in these 
nuclei, the FTBCSl and FTBCSl+SCQRPA predict similar results for the pairing gap and 
total energy. At the same time, for light systems, this deviation is stronger, therefore, the 
FTBCSl+SCQRPA offers a better approximation than the FTBCSl in the study of thermal 
pairing properties of these nuclei. 

The fact that the total energies and heat capacities obtained within the FT- 
BCSl+SCQRPA predictions agree reasonably well with the exact results for = 10 as 
well as those obtained within the finite-temperature quantum Monte Carlo method for ^^Fe 
shows that the FTBCSl+SCQRPA can be applied in further study of thermal properties of 
finite systems such as nuclei, where pairing plays an important role. Compared to existing 
methods, the merit of the present approach lies in its fully microscopic derivation and sim- 
plicity when it is applied to heavy nuclei with strong pairing, where the effect of coupling 
to SCQRPA is negligible so that the solution of the SCQRPA can be avoided. In this case, 
thermal pairing can be determined solely by solving the FTBCSl gap equation, which is 
technically as simple as the FTBCS one, whereas the exact diagonalization is impracticable 
(at T ^ 0). 

As the next step in improving the developed approach, we will include the effect of angular 
momentum in this approach. This study is now underway and the results will be reported 



in a forthcoming article 



Ml 



27 



Acknowledgments 



The authors thank Vuong Kim Au of Texas A&M University for valuable assistance. 
NQH is a RIKEN Asian Program Associate. 

The numerical calculations were carried out using the Fortran imsl Library by Visual 
Numerics on the RIKEN Super Combined Cluster (RSCC) system. 

APPENDIX A: FACTORIZATION OF {A]Aj>) 

The factorization of the screening factor {A'jAj') is not unique as it can be carried out 
in at least two ways, which lead to different results. In the first way, one can perform the 
mean-field contraction by using the Wick's theorem (WT) to obtain 

{AlA,,)-^^ ^ . (Al) 



In the second way, one uses the Holstein-Primakoff 's (HP) boson representation 35 1 



^j = 6jyi--^, A, = ^l-j^b,, M, = 2b]b,, (A2) 

with boson operators fej and bj to obtain 

{A]Af)^^ - 5jj'Vtjnj{l - 2nj) . (A3) 

The lowest order of the HP boson representation implies that operators A^^ and Aj are ideal 
bosons foj and 6j, respectively, i.e. setting Vj = 1 in Eq. ([6]). It is in fact the well-known 
quasiboson approximation (QBA), which is widely used in the derivation of the QRPA 
equations. The QBA leads to 

^ l^n'm = ^n'^^n, . (A4) 

As for the screening factor (^^.4],), it vanishes in these approximations. 

Using these results, we obtain the same form of Eq. (^U^ for the pairing gap, except that 
now e'j = ej, and the level-dependent part 6Aj from Eq. (15UI1 becomes 

SAf^ ^ 2GujVjnj , (A5) 
6Af ^ 2Gu,v, -Al-^.H2nj~l)n,] ^^^^ 
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FIG. 6: (Color on line) Level-weighted gaps for = 10 with G = 0.4 MeV as predicted by the 
WT (dashed), HP (dash-dotted), and QBA (thin dotted) approximations in comparison with the 
FTBCS (thick dotted), FTBCSl (thin solid), and FTBCSl+SCQRPA (thick solid) results. 

5 Af ^ ^ 2Gu,v, ^^ (1-^^ -^^) , ( A7) 

J i — 2nj 

which correspond to the approximations using the Wick's theorem, HP representation, and 
the QBA, respectively. 

The level-weighted gaps A obtained for = 10 and G = 0.4 MeV within these approx- 
imations are compared with the FTBCS, FTBCSl and FTBCSl+SCQRPA resuhs in Fig. 
El At T < 1 MeV, all three approximations, WT, HP, and QBA, predict the gaps close to 
the FTBCS one, but collapse at different Te, namely T^^^ < T™^^ < T^^ < T^t. At 
T ~ 1.2 MeV the HP gap reappears and increases with T to reach the values comparable 
with those predicted by the FTBCSl and FTBCSl+SCQRPA at T > 2 MeV. From this 
comparison, one can see that the mean-field contraction (1A1I) for {AjAj') includes only a 
tiny fraction of the QNF because it produces a finite gap at T > J'Jtbcs ^ g^g MeV, but 
this gap collapses again at Tj^'^ ~ 0.6 MeV. The HP boson representation, on the other 
hand, is able to take into account the effect of QNF at hight T leading to a finite gap at T > 
1.38 MeV, but fails to account for this effect at intermediate temperatures 0.55 <T< 1.38 
MeV. The QBA produces essentially the same result as that of the conventional FTBCS at 
low T with a slightly lower critical temperature T^^^ ~ 0.43 MeV. However, it causes a 
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negative A at T > 1.9 MeV. 
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